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Abstract 

Trajectory scaling functions are the basic element in the study of 
chaotic dynamical systems, from which any long time average can be 
computed. It has never been extracted from an experimental time 
series the reason being its sensitivity to noise. It is shown, by numerical 
simulations, that the sensitivity of the scaling function is to drift in 
the control parameters, and not noise. It is also explained how naive 
averaging of the orbit points may lead to erroneous results. 

The experimental study of chaotic dynamical systems presents us with 
complicated geometrical objects — strange sets — that have to be sim- 
ply characterized to be compared with theoretical predictions. The strange 
attractors are reconstructed from experimental time series through a pro- 
cedure known as phase space reconstruction [|, which determines the 
strange attr actor up to coordinate transformations. Therefore any charac- 
terization of the strange attractor must be independent of the coordinates 
used. Many different functions and sets of numbers have been proposed to 
characterize strange attractors, such as f(a) spectrum of singularities || 
and fractal dimensions, but the only complete characterization is the one 
given by the scaling function H, defined later on. There have been few at- 
tempts to extract the scaling function from experimental data |5], [f| , due 
mainly to its sensitivity to noise in the system. In this paper I will make 
explicit the difficulties and analyze a proposed extraction method: that of 
averaging the behavior of the system in the reconstructed phase space ||. 
I will show that the averaging procedure, as proposed, is not an effective 
procedure to extract the scaling function from an experimental time series. 
The main difficulty is that although averaging does reduce noise, it does not 
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reduce the main source of error in extracting the scaling function which is 
the detuning of the external parameters from the ones where theory makes 
its predictions. 

The scaling function a(t) was introduced by Feigenbaum and gives 
the local contraction rate of an asymptotically long periodic orbit after 
transversing a fraction t of the orbit. From it all other quantities of physical 
relevance can be explicitly computed. Scaling functions should be contrasted 
to other quantities that are extracted from dynamical systems, such as gen- 
eralized dimensions and f(a) spectrum of singularities. Although these 
quantities are invariants of the dynamical system (they remain unchanged 
if coordinates are changed), it is not possible to use them to compute all 
physically observable quantities such as the average energy dissipated in a 
chaotic circuit or correlations in the time series. The proof that the scaling 
function can be used to compute all physical averages was given by Sullivan 
H and also by Feigenbaum ||. 

The approach to understanding the effects of noise and systematic errors 
will be through the numerical simulation of circle maps. I will review the 
sine circle map in section [l] and give a few of the definitions that will be 
used later on. The various type of errors that hinder the extraction of the 
scaling function from time series are discussed in section |2[ Systematic 
errors will also be discussed in that section, as they are the major source of 
error in extracting the scaling function. The details on how to compute the 
scaling function are discussed in section ||; in particular I will concentrate 
on how to extract an approximation to the scaling function for the golden 
mean rotation number. All these sections are preliminaries to the results 
discussed in section |], where the sine circle map with noise in the parameters 
is explained. The surprising result is that very large noise levels have little 
effect on the scaling function when compared to systematic errors. In that 
section I also discuss a non-ergodic behavior of circle maps that occurs while 
averaging. 

1 Circle maps 

Maps of the circle occur whenever two oscillators are nonlinearly coupled. 
In general the asymptotic behavior of the coupled oscillators can be well 
described by a map that gives the difference in phase between them. For a 
circle map there are two relevant parameters: one which controls the ratio 
of the frequencies between the oscillators when they are uncoupled (u in 
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equation (|l])), and the other which controls the amount of coupling between 
the oscillators (k in equation @)). An example of a circle map that arises 
from the study of Hamiltonian systems is the sine circle map: 

k 

Xi+i = Xi + oj - — sin(27rxj) . (1) 

This is a map from the circle (parameterized from to 1) to itself, that 
is, all iterations of the map are computed mod 1. This map models the 
phase difference between two coupled oscillators. The interesting property of 
coupled oscillators is that they can mode-lock — while one of the oscillators 
executes p cycles, the other goes through exactly q cycles. The fraction p/q 
is the rotation number of the map and it represents the average fraction 
of the full range of the map transversed by each iteration. In general the 
rotation number is defined as 

p = lim — , (2) 

n— »oo fi 

with x n computed without the mod 1 after each iteration. If the strength 
of the coupling is non-zero, then there is a connected region in parameter 
space (k, u) where the rotation number is constant: the Arnold tongue. 

As the strength of the coupling between the two oscillators increases, 
larger ranges of uj are part of a tongue. At k = 1 almost all values of uj 
belong to some tongue and the map is said to be on the critical line. If 
the rotation number p of the map is an irrational number then the orbit 
of the map will be chaotic due to an (instant) period doubling cascade at 
the critical line. This is only proven for a class of irrational numbers with 
a particular number-theoretic property: if we expand the rotation number 
into a continued fraction expansion, then the terms of the expansion will 
not grow faster than a given power. If 

p= — = (0,1,0,2,...) (3) 



a± + 



0,2 H 



is an irrational number, then for constants C and 6, the coefficients a n of 
the expansion are bounded 

a n < C6 n . (4) 

The simplest proof of this chaotic behavior is through a renormalization 
group construction [|l(]] which is simplest for the golden mean irrational 
p g = (1,1,1,...). In what follows I will concentrate on the golden mean 
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rotation number. At this rotation number, the behavior of the map can be 
approximated by considering a sequence of maps with rotation number given 
by the approximants Q n /Qn+\ of p g obtained by truncating its continued 
fraction expansion. One finds that Qq = 1, Q± = 2, and Q n +i = Q n + Q n -i 
(the Fibonacci numbers). 

2 Noise 

Noise in a dynamical system can be present in many forms: in observations, 
in the state, or in the dynamics. If the system evolves deterministically 
under a map F r depending on a parameter r, but the position (state) is 
not measured accurately, then there is observational noise. This correspond 
to having the dynamics = F r (xi), but observing X{ + £j rather than 
Xi, where £j is a random variable (noise). The system may also evolve 
stochastically. In this case the noise may change the state at each time step 

Xi+i = F r {xi) + & , (5) 

or it may change the dynamics 

x i+ i = F r+ £.(xi) . (6) 

Combinations of all three types may occur and in general all are present in 
a laboratory experiment. 

The scaling function is very sensitive to noise and to the parameter values 
of the map, which has made it difficult to extract it from experimental data 
or even from numerical simulations. As the outcome of most experiments 
with chaotic systems is a time series, I will concentrate on how scaling 
functions are extracted from them. The simplest method to eliminate the 
error in the time series is by averaging it over several periods. Even though 
averaging can diminish observational and state noise, it does not change 
the fact that there are drifts in the experiment that lead to systematic 
errors. As we will see, averaging over periods does little to diminish the 
error in computing the scaling function, as it does not change the errors 
made in tuning the parameters to the golden mean. Because observation 
noise and state noise can be made small in an experimental setup (by care 
in the experiment or by period averaging), I will only consider the effects of 
dynamical noise and systematic errors. 

In experiments with systems at the borderline of chaos, systematic errors 
are the largest. When collecting the data for the circle map the experimen- 
talist has to tune the parameters so that the rotation number is exactly the 
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golden mean. The golden mean is an irrational number and its associated 
tongue has no width, which makes the tuning only approximate. Then to 
extract the scaling function, or even a simpler thermodynamic average such 
as the f(a) spectrum of singularities, the longest possible data set must be 
collected, which implies that the parameters must be kept at the golden 
mean for a long time. The time scale is determined with respect to the 
natural frequency of the system, if it is a self-oscillator, or with respect to 
the external frequency, if it is a forced oscillator. In a typical experimental 
setup the parameters cannot be kept tuned to golden mean rotation number 
and a slow drift in the parameters can be detected. This drift is interpreted 
in the experiment as a systematic error. 

If the rotation number is determined from the Fourier spectrum, then 
its accuracy is low. If N data points are used in computing the spectrum, 
then the rotation number, which is a frequency, is known to an accuracy of 
order 1/JV. Better techniques for computing the rotation number have been 
developed which take into account that the orbit points have a well defined 
ordering around the circle. With these techniques it is possible to determine 
the rotation number to an accuracy of order 1/N 2 

To convey an intuition on how sensitive an experiment can be to drift 
consider the Rayleigh-Bernard convection experiment performed in a mix- 
ture of 3 He and 4 He at mili-Kelvin temperatures [12, |13|| . The data from 
this experiment was collected for several days without interruption and the 
relevant control parameter — temperature — was kept tuned to the golden 
mean value to within 1 part in 10 5 . Nevertheless 24 hour fluctuations on 
the rotation number could be seen while the laboratory air-conditioner was 
turned off. Once the laboratory temperature was regulated other fluctu- 
ations on the scale of an hour could be detected with the 1/N 2 method. 
Variation of the position in parameter space seems unavoidable in any ex- 
perimental setup. 



3 Scaling Function 

The scaling function for a circle map is computed at a quadratic irrational 
which has a periodic continued fraction expansion. As the inflection point, 
xo, is iterated it rotates on average p, and the points xo,xi,. . . of the orbit 
delimit a series of intervals or segments along the circle. The endpoints of 
these segments are not two successive iteration points, such as Xt and xt+i, 
but depend on how many times the initial point has been iterated (see figure 
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Figure 1: The segments used for the construction of the scaling function are 
determined from the first iterates of the map. The numbers n on the top 
segment label the iterates x n of the map. By considering the orbit points 
separated by a Fibonacci number, the segments can be arranged in levels. 



|H). If the number of iterations is a Fibonacci number, then the orbit points 
can be arranged in groups (or levels) that recursively subdivide the circle 
into smaller and smaller segments. An example of the subdivision processes 
is given in figure ffl. The first 13 orbit points of the golden mean trajectory 
are indicated in the figure. Notice that the segment A is delimited by the 
orbit points xo and xq„, and that its location alternates to the right and left 
of the point xq. The other segments of the level are determined by mapping 
the segment A around the circle. For universality, this construction is not 
to be carried out with the actual map, but rather with a Q n iterate of the 
map 0, m. 



For the case of a simple repeating number in the continued fraction 
expansion, such as the golden mean, the segments are given by 

A^ = \x s -x Qn+s \ (7) 

from which we can define the values assumed by the scaling function at 
different points 

a ( n ) = (8) 

A^[s<Q n }+A%\s>Q n ] 



6 



A square brackets [16] evaluates to one if the expression within them is 
true and zero otherwise, so that the denominator of the expression chooses 
one of the segments, or A^_Q n , as appropriate for the segment on 

the numerator (see figure [l]). An approximation to the scaling function 
is obtained by the concatenation of Q n short steps of length l/Q n and 

(n) 

height as in ascending order of s. This defines a function from the unit 
interval to itself. The approximation in terms of steps of constant height 
is a reasonable approximation because the variation in height of the steps 
diminishes exponentially fast as the number of the level n increases. The 
construction of the continuous (and also differentiable) almost everywhere 
scaling function is 

where [x\ is the function that gives the largest integer smaller than x. When 
evaluating the scaling function from a map with the parameters different 
from the golden mean rotation number, then there is another limit involved: 
that of approaching the irrational number winding number. The two limits 
do not commute, and the irrational winding number must be approached 
before the limit to an infinite number of levels. In practice the irrational is 
approached as well as possible. Also, for universality, the scaling function 
must be computed in a neighborhood of the inflection point. In practice 
the problem is bypassed by taking Qq to be not 1, but a larger Fibonacci 
number (see reference [11]). This follows from the properties of circle maps 
with a golden mean rotation number. Under composition by a Fibonacci 
number of times the orbit points accumulate around the starting point for 
the iterations, which is taken to be the inflection point. 

The limit @ to compute the scaling function has to be approximated 
in practice with a large enough n. From experiments it has proven feasible 
to extract the scaling function which has 5 steps. This approximation is 
plotted in figure [2] together with the limiting scaling function. 



4 Simulations 

To understand the sources of error in determining the scaling function, I will 
compute it from an orbit of a map that is not exactly at the golden mean 
rotation number but at one of its continued fraction approximants with a 
length typical of what is obtained in laboratory experiments. In particular I 
will consider the orbit with rotation number 21/34, which in the circle map 
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Figure 2: Three five step approximation (in black) to the limiting scaling 
function (in gray). 

occurs at u = 0.606439 on the critical line (k = 1). To simulate the effects 
of noise fluctuations the control parameters k and u of the sine circle map 
will be slightly varied at each iteration. Both k and u will be replaced by 



where Ak and Aco are the strength of the fluctuations and and s, are 
random numbers uniformly distributed in the interval from —1 to 1. Notice 
that k and u remain fixed during the random process, and for the uniform 
distribution, represent the average values of ki and Wj. An orbit from the 
noisy sine circle map is generated by 



If the average parameter values are within the 21/34 tongue, then the map 
is iterated a few hundred times before any orbit points are used to compute 
the scaling function. If the average parameters are not within the tongue, 
then the map is started at the inflection point. 

The orbit is averaged according to the procedure proposed by Belmonte 
et al. ||, where points that are nearby in coordinate space are averaged 
together and coalesced into a single orbit point of an averaged periodic 



ki = k + riAk 
oji = to + SiAui 



(10) 



x i+ i = Xi + uji- - 1 sin(27ra;i) . 
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orbit. If the average parameter values of the map are within the tongue of 
the rotation number being considered, the group of points to be averaged 
can be unambiguously distinguished for errors as large as A/c = Au> = 0.05, 
an error much larger than in most experiments. If the average parameter 
values are outside the tongue, then the number of groups to be averaged 
will depend on the length of the data set. 

The first observation from the numerical simulations is that small errors 
can lead to largely distorted scaling functions. In figure |3|a the scaling func- 
tion for a short orbit (21/34) at parameters k = 0.9999 and to = 0.6063, 
which is close to the superstable point of the tongue, is compared with 
the theoretical curve. The amplitude of the error fluctuations are small 
(Ak = Auj = 10~ 4 ) which keeps the map parameters within the tongue. 
In this case there are large deviation from the theoretical curve. In figure 
||b, for the same orbit length, the scaling function is computed with fluctu- 
ation noise 100 times larger, but with parameters {k = 1.0 and uj = 0.6066) 
closer to the golden mean critical point. The difference between the scaling 
function obtained from the short orbit and the theoretical curve is smaller 
than in figure ||a. This at first seems paradoxical: the curve with larger 
fluctuations is closer to the theoretical curve than the curve with smaller 
fluctuations. 

To quantify the differences between the theoretical and short period 
scaling function the L 1 norm can be used. This norm is proportional to the 
area in between both curves. If a{t) is the theoretical scaling function with 
five steps, and a (t) is the scaling function obtained from the orbit with 34 
points, also with five steps, then the error between them is defined as 

1 f 1 

e(a, a ) = — dt \a(t) - a (t)\ , (12) 

Co JO 

where Co is a normalization constant. The constant is chosen so that error 
between the theoretical scaling function and the one obtained from a short 
orbit at the irrational winding number is one. The constant cq is computed 
to be 0.01121. With this norm the error between the curves in figure |3|a is 
2.26 and between the the curves in figure |3|b is 0.53. 

The systematic error constitute the larger source of error. This can be 
verified by plotting the error between the scaling function obtained at a 
point away from the golden mean rotation number and another point along 
the critical line. The further the rotation number is from the golden mean, 
the larger the difference between the two scaling functions. In figure || the 
inflection point is iterated for 34 times; from this orbit a five step scaling 
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Figure 3: Different type of errors lead to different scaling functions. In both 
figures the theoretical curve is indicated by a solid line. In figure (a) the 
fluctuation errors are small, but the scaling function deviates largely from 
the theoretical curve. In figure (b) the fluctuation errors are large, but the 
scaling function deviates only slightly from the theoretical curve. 



function is computed, which is then compared to the asymptotic five step 
scaling function. In the figure the rotation number is measured from its 
departure from the golden mean rotation number in units of the width of 
the 21/34 tongue. In actual units of the map the horizontal axis ranges 
from 0.60638 to 0.60685, which is four times the width of the 21/34 tongue. 
According to the plot, the error is smallest when the rotation number is 
closest to the golden mean, and increases as one departs from it on either 
side. The exact minimum in the error curve does not coincide with the 
golden mean because of finite size effects in computing the scaling function. 
For the error curve to be a smooth function of the rotation number it is 
necessary that all orbits start at the same point, the inflection point in this 
case. 

The plot of figure || was obtained from iterating a map without fluctua- 
tions in the control parameters. One would expect the the results obtained 
without the error are those that would be obtained had the map with error 
been iterated and averaged a large number of times. This is the case if the 
points are averaged according to their time index, that is, for an orbit of 
period P, the average over the fluctuations of the k-th point of a periodic 
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Figure 4: The error in approximating the scaling function by an orbit gener- 
ated from a map with systematic error. The winding number uj is measured 
in units of the width of the 21/34 tongue away from the golden mean and 
the error is the area in between the curves. The arrow shows the location 
of the point for the 21/34 cycle 

orbit are computed from 

(x k ) = lim - V Xk+iP ■ (13) 

n^oo n * — ' 
l<i<n 

But this may not be the average that is computed in an experiment. Some- 
times it is simpler, or consistent with time delay coordinates, to average 
the points that are close to each other in time delay space (this was the 
procedure adopted in reference ||). In table [l] an orbit for a map at the 
superstable point has been iterated with a small error (Ak = Aw = 10~ 3 ). 
The map is iterated while the parameters fluctuate. Each point is compared 
with the exact orbit and iterated points that come close to the same exact 
orbit point are averaged together. From the averaged orbit the five level ap- 
proximation to the scaling function is computed and used to determine the 
error associated with the orbit by comparing it to the scaling function with- 
out noise. By "without noise" , I mean the scaling function that is obtained 
by iterating the map with the average parameters values of the simulation 
with noise. The table shows the results of longer and longer averages. At 
first the error diminishes, but as the number of samples increases the error 



11 



samples 



error 



10° 
10 1 
10 2 
10 3 
10 4 
10 5 



1.99065 
0.21849 
0.53230 
0.49462 
0.52383 
0.51801 



Table 1: Error between the scaling function computed with and without 
noise. The noisy map has fluctuating parameters with average at the su- 
perstate point of the 21/34 tongue. The averaging is done in coordinate 
space. As the number of samples increases the error does not go to zero, as 
would be expected. 

appears to remain constant. The conclusion from the table is that one has 
to be careful that the limits involved in the averaging procedure are well 
defined and converge to ones expectations. 

5 Conclusions 

From the numerical simulations one sees that even large errors can have 
little effect on the extraction of the universal scaling function, provided that 
the parameters of the system are well tuned to the golden mean rotation 
number at the transition to chaos, as can be seen from figure ||. The error 
in computing the scaling function depends on how close the parameters are 
to the transition point to chaos, a quantity that is difficult to control in 
experiments, as they are invariably subject to drift. The drift comes from 
the conflicting requirement of tuning the parameters to the smallest possible 
tongue (and therefore at the limit of instrumentation) and of obtaining the 
longest possible times series. 

Also from the numerical simulations the perils of averaging the orbit in 
coordinate space where pointed out. The noise in the system causes the 
orbit to land close to the "wrong" group of points for its phase within the 
period, which leads to the non-convergence of the averaging procedure. I 
have no mathematical proof for this lack of convergence, but table |] gives 
numerical evidence towards the result. 
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It hardly seems worthwhile to extract the scaling function given all these 
difficulties, specially given that the f{a) spectrum of singularities seems very 
robust to noise and simple to extract from experimental time series. It also 
appears to give an infinity of scales for the problem, just as the scaling 
function. The difficulty with this argument lies with the error bars of the 
spectrum of singularities. With error bars of the order of 1%, the spectrum 
of singularity is equivalent to just three of the values of the scaling function 
[17 1; the spectrum of singularities does not give individual scales but mixes 
them all into one function. In order to extract further information from the 
spectrum of singularities the errors would have to be reduced well below 
the 1% level, which does not seem possible even in numerical simulations. 
Contrast this with the scaling function. There every different scaling (the 
ai™^) can be individually and independently extracted. 

Scaling functions are the fundamental objects in the study of low dimen- 
sional dynamical systems, and it is surprising how little attention they have 
received in the literature, both theoretically and experimentally. If they are 
to be extracted from experimental time series new techniques will have to 
be developed to control the systematic errors. 
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